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(57) Abstract: A system and method of creating a filter for 
use with locally dense seismic data is disclosed. The method 
includes obtaining survey geometry characteristics from a lo- 
cally dense seismic survey A filter is designed which uses 
spatial derivatives of the wavefield of order between (1) and 
the maximum order of spatial derivatives of the wavefield that 
can be estimated within a group. The filter can be designed so 
as to separate up/down going components, p/s components, 
or both up/down and p/s components. Partial derivatives in 
space and time of the wavefield can be calculated, using, for 
example, a laylor scries expansion as an approximation. The 
seismic data is filtered by combining estimated near surface 
material properties, the seismic data, and the calculated par- 
tial derivatives. 
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System and Method for Seismic Wavefield Separation 

FIELD OF THE INVENTION : 

The present invention relates to the field of 
5 seismic data acquisition and processing. In 

particular, the invention relates to a system and 
method for seismic wavefield separation. 

BACKGROUND OF THE INVENTION : 

10 Optimal processing, analysis and 

interpretation of land seismic data ideally require 
full information about the wavefield so that the 
wavefield can be separated into its up- and down-going 
and P- and S-components as well as determining phase 

15 and polarity. For 3C acquisition of land surface 

seismic data it is common practice simply to interpret 
the vertical component as the P- section and the 
horizontal components as SV- and SH-sections. This 
"traditional" P/S interpretation is exact for vertical 

20 arrivals. However, as energy is incident away from 
normal incidence angles, this approximation breaks 
down, both because of projections on to all components, 
but also because of a non-unity reflection coefficients 
and mode-conversions at the free-surface. 

25 Exact analytical filter expressions for 

wavefield separation have previously been derived by 
for instance Dankbaar, J. W. M. , 1985, Separation of P- 
and S-waves: Geophys. Prosp., 33, 970 - 986, and these 
have been applied to seismic data in conventional 

3 0 recording geometries. Unfortunately, statics problems 

severely limit the practical use of these wave-equation 
based techniques . 
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SUMMARY OF THE I3WENTION : 

Thus, it is an object of the present 
invention to provide a filtering technique that lends 
itself to be applied within local densely deployed 
5 single sensor receiver groups where statics are 
substantially constant. 

It is a further object of the present 
invention to provide a filtering technique that can be 
implemented efficiently directly in the spatial domain. 

10 According to the invention a method of 

creating a filter for use with locally dense seismic 
data is provided. The method includes obtaining survey 
geometry characteristics from a locally dense seismic 
survey designed to record characteristics of an elastic 

15 or acoustic wavefield. The seismic survey is made up 
of a number of groups of receivers, with each group 
comprising at least three receivers densely spaced from 
each other. According to the method, a filter is 
designed which uses spatial derivatives of the 

20 wavefield of order between 1 and the maximum order of 
spatial derivatives of the wavefield that can be 
estimated within a group. The filter is designed such 
that when combined with data from within a single 
group, the filter separates components of some or all 

25 of the wavefield arriving at the single group. 

The filter can be designed so as to separate 
up/down going components, p/s components, or both 
up/down and p/s components. 

The densely spaced receivers in the group are 

3 0 preferably spaced apart such that statics in the 

portion of the wavefield of interest are substantially 
constant. More preferably, each of the densely spaced 
receivers in the group are spaced apart by about 2 
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meters or less, or by a distance of about one fifth the 
shortest wavelength of interest or less. 

Partial derivatives of the wavefield are also 
preferably calculated, and this can be done using a 
5 taylor series expansion as an approximation. According 
to the invention, the seismic data is preferably 
filtered by combining estimated near surface material 
properties, the seismic data, and the calculated 
partial derivatives (both in space and time) . 

10 The filter can also be used to separate 

surface waves or airwave induced ground motion from the 
seismic data. The free surface condition can be used 
to convert vertical derivatives of the wavefield to 
horizontal derivatives of the wavefield. 

15 According to the invention, the seismic 

survey is performed primarily for the purpose 
hydrocarbon reservoir exploration, evaluation or 
characterisation, although other uses can be made. 

The invention is applicable where the near 

20 surface velocities are isotropic or anisotropic. 



BRIEF DESCRIPTION OF THE DRAWINGS : 

Figures la-c show various illustrations of 
25 the elastodynamic representation theorem, according to 

the invention; 

Figure 2 shows a plot of the normalised real 

part of the filter for removing free surface effects 

from divergence in equation (36) ; 
30 Figures 3a and 3b show an example of a 

seismic survey having locally dense receiver groups, 

according to the invention; 
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Figures 4a and 4b show the divergence as 
calculated directly by numerical methods; 

Figures 5a-d illustrate sections obtained 
from v 3 measurements only; 

Figures 6a-d show the results of using a 
wavefield separation technique according to a preferred 
embodiment of the invention; 

Figure 7 is a schematic illustration of a 
seismic data acquisition and processing system, 
according a preferred embodiment of the invention; and 

Figure 8 shows steps involved in the 
filtering technique, according to preferred embodiments 
of the invention. 

DETAILED DESCRIPTION OF THE INVENTION : 

According to the invention a new approach is 
provided for P/S separation of land surface seismic 
data and for removing the effects of the Earth's free 
surface (i.e. up/down separation). By converting 
vertical spatial derivatives to horizontal derivatives 
using the free-surface condition, the methodology can 
make use of locally dense measurements of the wavefield 
at the free surface to calculate all spatial 
derivatives of the wavefield. These can in turn be 
used to compute divergence (P-waves) and curl (S-waves) 
of the wavefield at the free surface. 



preferably removed through an up /down separation step 
using the elastodynamic representation theorem. This 
results in expressions where spatial filters are 
convolved with recorded data. The filters can be 
successfully approximated so that they fit the local 



The effects of the free surface are 
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dense acquisition pattern used for the P/S separation 
step. In particular, the simplest approximate 
expression for up-going P-waves consists of two terms. 
The first term preferably corresponds to the divergence 
5 in the presence of the free surface scaled by a 

material constant. The second term preferably is a 
time derivative of the recorded vertical component 
scaled by a material constant. Hence, a correction is 
added to the u traditional" P- interpretation through the 
10 first term, which improves accuracy for incidence 
angles outside normal incidence. 



Method for Estimating Seismic Material Properties" (UK 
Patent Application No. 0003410.8) filed on 15 February 

15 2000, incorporated herein by reference, a method is 

provided than can use the volumetric recordings of the 
wavefield to invert for P and S velocities in the Earth 
in the neighbourhood of a small, closely-spaced array 
of receivers . "Volumetric recording" refers to an 

20 array that approximately encloses a volume of the 
Earth. The quantities estimated are the effective 
velocities of the P- and S-components of the wavefield 
at any point in time. Hence, these will vary with both 
wave type and wavelength. If estimated for the near 

25 surface Earth structure, such velocities may be useful 
for statics estimation, or for separation of the 
wavefield into up- and down-going components as 
described herein. 



30 for land seismic recording will be discussed first with 
reference to the fully anisotropic case. The elastic 
constitutive relation relates components of the stress 



In UK Patent Application entitled "System and 



Implications of the free surface condition 



WO 01/53854 



- 6 - 



PCT/GB01/00186 



tensor ay in a source free region to the strain tensor 
components €y 

Qij = Cijkjejci , (1) 

5 

where cga are the elastic stiffnesses. Index 
values 1 and 2 correspond to horizontal coordinates xj 
and X2 whereas index value 3 corresponds to the 
vertical direction downwards xs . Using the Voigt 
10 notation, equation (1) can be written as: 
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where C is the symmetric stiffness matrix 
15 with 21 independent components: 
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(3) 



The strain is related to particle velocity 

2 0 V,- through 



(4) 
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where 9i denotes spatial derivatives in the 
xi, X2 or X3 directions, and the dot denotes a time 
derivative. In land acquisition; spatial distributions 
5 of 3C receivers at the surface allow us to compute 
horizontal spatial derivatives of particle velocities 
(or time-derivatives thereof if particle acceleration 
is recorded, etc.). The only information missing in 
order to know the wavefield completely is the vertical 
10 derivatives of the recorded wavefield. 

The free-surface condition at the Earth's 
surface gives us three additional constraints: 



= 0 , (5) 

15 

which are sufficient to allow us to compute 
the remaining velocities given that we make some 
independent estimate of the relevant elastic 
stiffnesses . 

20 In addition, constraints on the relation 

between individual elements of the stress and strain 
tensors could possibly be used to correct for coupling 
or to compute near-surface properties. 

The case of isotropic material properties in 

25 the near-surface environment is of particular interest 
in the land surface seismic case. Using the Voigt 
notation, the stiffness matrix takes the following 
form: 
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c = 
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(6) 



where X and ju are the Lame constants. 
The constraints imposed by the free-surface 
5 condition (5) become: 



A + 2pL 

d 3 t> 3 , 



(8) 
(9) 



where the horizontal derivatives on the right 
10 hand side are known from the surface measurements. 
Note that the material properties only occur in 
equation (7) and not in equations (8) and (9) . 

Divergence and curl of a wavefield at a free 
surface overlaying a homogeneous isotropic half space 
15 will now be discussed. This discussion is in the 
context of an isotropic media. The elastic wave 

equation for particle displacement u can be written as: 

i 

pu = f + (A + 2p)V(V ■ u) - fiV X (V X u) , (10) 

20 

where p is the density and f denotes a 
distribution of body forces. Lame's theorem states 
that there exist potentials <I> and W of u with the 
following properties: 



WO 01/53854 

# 



- 9 - 



PCT/GB01/00186 



u = V$ + V x * , (11) 

V* = 0, (12) 

$ = i +c 2v 2 ^, (131 
P 

4> = ^ + 4v 2 * 3 (i4) 



where c a and cp are the P- and S-velocities , and (j> and y/ 
5 are potential components associated with the body 

force. An elastic wavefield u can thus be decomposed 
into its P- and S-wave components, V<2> and VxT, 
respectively. Equations (11) and (12) yield: 

V-u = V 2 $, (15) 
io Vxu = VxVx* (16) 

By measuring the curl and the divergence of 
an elastic wavefield we can thus measure the P- and S- 
wave components separately. 

15 An acquisition pattern comprising tetrahedra 

of 3C measurments can be used to achieve the separation 
of the wavefield into its curl- and divergence-free 
components. See, e.g. Robertsson, J.O.A., and Muyzert, 
E., 1999, Wavefield separation using a volume 

20 distribution of three component recordings: Geophys . 

Res. Lett., vol. 26, 2821-2824, incorporated herein by 
reference. Such acquisition patterns are described in 
further detail in UK Patent Application No. 9921816.6, 
incorporated herein by reference. All spatial 

25 derivatives of the wavefield components can be 
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calculated. Therefore, divergence and curl can be 
calculated explicitly from surface measurements only. 
Equations (7) , (8) and (9) give us the following 
expressions for divergence and curl of particle 



5 velocity at a free surface: 

(v " v) = jf^^ 1 * 1 * 92 ^ 1 (17) 

(Vxv)i = 23^3, (18) 

(Vxv) 2 = -26 1 v 3 , (19) 

(Vx v) 3 = dit? 2 - d 2Ul . (20) 

First we note that the ratio £ — = 2{cp/c a ) 



(may be frequency dependent) scales the expression for 

10 the divergence in equation (17) . 

Second, equations (17), (18), (19) and (20) 
contain both the up-going and down-going parts of the 
wavefield. This includes mode conversions at the free 
surface. For instance, the divergence given by (17) 

15 contains not only the desired up-going P-waves, but 
also, the down-going P-to-P reflection, down-going S- 
to-P conversions, etc. Moreover, a plane P-wave which 
is vertically incident on the free-surface will have 
zero divergence (up- and down-going parts interfere 

20 destructively) . Removing the effects of the free 

surface is therefore regarded as an important step in 
the preferred P/S separation technique. 

A technique for Up/down wavefield separation 
using the elastodynamic representation theorem will now 

25 be discussed. The elastodynamic representation theorem 
or Betti's relation can be derived from the equation of 
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motion and the elastic constitutive relations using 
Gauss' theorem. Suppose that we have a volume V 
enclosed by a surface S, and that we wish to calculate 
the displacement of a wavefield u at a point x' in V. 
5 The displacement u(x') is directly related to the stress 
o and displacement along S f sources in V, as well as 
the displacement Green's tensor Gy and the stress 
Green's tensor 2p functions between S or source points 
and x' : 

^(x ; ) - £(t<x)G^ (21) 
10 Jv tw 



where n is the normal unit vector to S, U(x) 
is the traction across 5 at point x and fi is the i th 
component of the force from sources within V. 
15 Equation (21) shows the frequency domain expression of 
the representation theorem. In the time domain we 
would also have convolutions in time. The traction t 
along S is defined as 

20 t=n <7. (22) 

Now, suppose that the volume V consists of a 
homogeneous elastic medium. Figures la-c show various 
illustrations of the elastodynamic representation 
25 theorem, according to the invention. First, consider 
the closed surface Si+S 2 illustrated in Figure la. In 
this case we assume that we have a source located 
outside V so that it is recorded along Si as an up- 
going mode. Sommerf eldt 1 s radiation condition implies 
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that the contribution from S2 vanishes at infinity. 
Hence, the closed surface integral in equation (21) can 
be replaced by the open surface integral: 

u„(x') = / {t i{x)G in {x, x') - ti;(x)hjZ jin (x, x')) dS(x) . (23) 



The expression in equation (23) represents an 
up-going wavefield at x' since the field is up-going at 
10 Si. 

Next consider the situation depicted in 
Figure lb. If we in a similar way put the source above 
Si outside V, we obtain an analogous expression to 
equation (23) for the down-going wavefield. 

15 Finally, we consider the situation shown in 

Figure lc. This time we assume that' Si follows the 
Earth 's surface topography, that no sources (primary or 
secondary) are located in V, and that x' is located 
inf initesimally below the land surface. Again, the 

20 left and right edges of the surface S are at infinity 
and do not contribute to the integral . The 
representation theorem therefore becomes: 



In equation (24) we have used particle 
velocities v t instead of displacements u t since this will 
be more convenient to work with. From the above 
30 derivation it is now clear that the integral 

contribution over Si corresponds to the down-going 
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field at x' . Moreover, for a horizontal free surface 
we know that 0/3=0 (equation (5)). Since the total 
field consists of a sum of up-going and down-going 
waves, the up-going particle velocity field just below 
5 the surface is therefore given by 

v n (x ( ) = v n (x') + / ^(x)S 3) „(x J x / )(i5 , (x) . (25) 

Hereinafter the tilde denotes a wavefield 
10 with up-going waves only, i.e. having removed the 
effects of the free surface from the wavefield. 

According to a preferred embodiment of the 
invention, the elastodynamic representation theorem can 
be used to extract the desired up-going wavefield from 
15 surface seismic recordings. For a discussion in the 
case of seismic data from ocean bottom cables, see 
Published UK Patent Application GB 2 333 364 A, 
incorporated herein by reference. 

In the following discussion, co is the angular 
20 frequency, k a = o)/c a and kp = (o/cp are the P- and S- 

wavenumbers, c a and cp are the (isotropic) P- and S- 
velocities, and p is the density. 

The elastic displacement Green's tensor in an 
isotropic homogeneous medium is: 

25 

Gi„(x 9 id s v) = ^(^^n^(x,x^c;)+9^5 r (^(x 3 x^ v)-g a (x, x'»)) , (26) 



30 



Where g a and gp are the P- and S -Green's 
functions, and S in is Kronecker's delta. For 
homogeneous media in 3-D these are given by: 
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The corresponding isotropic stress Green's 
5 tensor is given in terms of the displacement Green's 
tensor as : 

= XSijd k G kn + ^{dfiin + 3,"Gy„) . (28) 

10 A flat horizontal free surface is assumed 

here. Expressions are therefore needed for Green's 
functions such that x = (^,X3) and x' =(£',*' 3) are located 
very close to each other in the ^-direction, whereas 
the horizontal separation can be arbitrary. Hence, X3- 

15 x'z- e -> 0 + and the Green's tensors are shift invariant 
with respect to the x\ and Xi directions. 

Since the "free-space" Green's tensors Gu and 
liijk in equations (26) and (28) are shift invariant, the 
expressions for the up-going wavefield given by the 

20 representation theorem (25) can be written as spatial 
convolutions denoted by * : 



*i« = |(t*(x) t*(x)), (29) 
*M = (30) 
*>3« = ^Wx) + ^(x)* Vl (x) + ^(x)*t; 2 (x)). (31) 



25 The filters in equations (29), (30), and (31) 

are: 
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F% = 2S 33 ,(s 3 = xj) = - 2d v ((l + 2k->d&)g a + 2£ 2 afo) , (32) 
?* = 2S, P a(jC3 = ij) =2^ 2 a,(-2^ ft +(^+2^)^), (33) 

where the subscripts v and £ denote either 
index 1 or 2 corresponding to horizontal coordinates X\ 
5 and Jt2- 

Using the Green's function in equation (27), 
we can obtain explicit expressions for the filters in 
equations (32) and (33) that can be implemented 
directly. In the ^fc-domain these are: 

0 

F+ = ^(l-a^'^ + W)) ■ 
= 1M i 1 - 2k p' + *W) ■ 



(34) 
(35) 



where kl a) = ^jk^-k^k^ is the vertical 

wavenumber for P-waves and fc^ = ^jk^ -k^k^ is the 

15 vertical wavenumber for S-waves. 

A preferred technique for P/S separation of 
surface seismic data will now be discussed. We can 
calculate the divergence and curl of the up-going waves 
by taking spatial derivatives of equations (29), (30) 

20 and (31) . Some care must be taken here since this 
involves spatial derivatives of the stress Green's 
tensor E^. In this fashion, the following expressions 
are obtained: 
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(V-v) = ^ lUl + a 2 v 2 ) + ^(x)*« 3 ( X ) 1 (36) 

(Vxv) a = ^T7 3+ ^ i Vx ^(y)*t ?1 (x)+^^( 5 c)*« J (x) J (37) 
(Vxv) a = -^ U3 +^ x ^(x)* yi (x)+^ 4 v>: ^(x)*t; a (x) ) (38) 
(Vx v) 3 = fava - dith) . (39) 

Again, using the Green's function in equation 
(27) , we can obtain explicit expressions for the 
5 filters in equations (36), (37) and (38). In the fk- 
domain these are: 



(40) 



r 3 



1>2 



(42) 



(43) 



10 The traditional way for P-wave interpretation 

is to simply look at the recorded V3 component . This is 
exact for vertically propagating P-waves . As we saw 
above, the divergence on the other hand is zero at the 
free surface. For horizontally propagating P-waves, the 

15 situation is reversed; V3 is zero whereas the 

divergence exactly contains the P-waves. From equations 
(36) and (40) we see that the correct expression 
combines v 3 with the expression for divergence in the 
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presence of the free surface so that it holds true for 
all incidence angles . 

Equations (40) to (44) are all functions of 

the material properties c a and cp at the receiver 
5 location. Estimation of these parameters using a full 
wave equation approach with dense surface measurements 
consistent with that taken here are described in UK 
Patent Application entitled "System and Method for 
Estimating Seismic Material Properties" (UK Patent 

10 Application No. 0003410.8) filed on 15 February 2000. 

An example of a preferred implementation of 
wavefield separation filters will now be discussed. The 
filters derived above can in theory be implemented 
directly which yield expressions that are exact for 

15 homogeneous media with a flat surface. These 

expressions would remove all down-going and evanescent 
wave types including ground-roll. 

However, the filters are slowly decaying and 
contain some complicating factors. High-order factors 

20 of ki and k 2 correspond to high-order derivatives in 

the spatial domain. The most severe problems arise from 
the factors k 3 (a> and k 3 (P> . These terms do not correspond 
to any straightforward implementation in the spatial 
domain. When they occur in the denominator they 

25 introduce a pole in k a and kp respectively. 

One straightforward filter approximation is 
to make Taylor expansions around = 0 in the 
wavenumber domain (factors of k^ correspond to spatial 
derivatives) . The lowest order terms in the Taylor 

30 approximations to equations (40), (41), (42), (43) and 
(44) are: 
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t£" } « - i k -fhhW - K') , (45) 

P (7XV) S „ i^__i. (2fc , + fccAt)i (4g) 

fij" 1 ' w -i^-- (49) 



Figure 2, shows a plot of the normalised real 
part of the filter for removing free surface effects 
5 from divergence in equation (40) . The solid line shows 
the exact filter (equation (40)). For wavenumbers 
beyond the pole corresponding to horizontally 
propagating P-waves, the filter is imaginary 
(evanescent modes) - Figure 2 also shows a plot of the 
10 zeroth-order (dashed line) , and the first-order (dash- 
dotted line) Taylor approximations as given by equation 
(45) . 

The application of a preferred wavefield 
separation approach to synthetic data will now be 

15 discussed. A reflectivity code was used to test the 

wavefield separation approach according to a preferred 
embodiment of the invention. See e.g., Kennett, B. L . , 
1983, Seismic wave propagation in stratified media: 
Cambridge University Press, Cambridge. This was chosen 

20 as opposed to for instance finite differences since up- 
and down-going wavefields can be calculated separately. 
Moreover, the quantities right below an interface can 
be obtained accurately. The output from the 
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reflectivity code are particle velocities and 
divergence of particle displacement. 

Figures 3a and 3b show an example of a 
seismic survey having locally dense receiver groups, 
5 according to the invention. Figure 3a is also the 
model used to test the wavefield separation approach 
according to a preferred embodiment of the invention. 
In the model, a point-source was used to generate 3-D 
synthetics. The source consisted of a 50 Hz Ricker 
10 wavelet and was of explosive type (radiating P-waves 
only) . 

As shown in Figure 3a source 202 is located 
100 m below earth surface 204. The wavefield was 
recorded at 25 m inter-group spacing from 0 m offset to 

15 450 m. Geophone groups 206, 208 and 210 are shown as 
the first three groups in a series of groups, each 
spaced apart by 25 meters. At each recording location, 
each receiver group comprises four 3C geophones were 
spaced evenly at 0.5 m in both horizontal directions. 

20 Also shown in Figure 3a is a reflective subterranean 
surface 212 that is located 150 meters below earth 
surface 204. 

Figure 3b is a plan view showing an example 
of geophone layout within a single receiver group 210. 

25 Receiver group 210 comprises four geophones 250, 252, 
254, and 256 located on the ground surface. Receivers 
250 and 254 are spaced apart about 0.5 meters from each 
other, and receivers 252 and 256 are spaced apart about 
0.5 meters from each other. Also shown in Figure 3b is 

30 a dashed line 260 on which the other receiver groups 
are positioned, for example, receiver groups 206 and 
208. 
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Advantages of moving to smaller inter- 
receiver spacings include increased accuracy in the 
estimated derivatives, and greater validity in the 
assumption that the material properties are the same in 
5 the vicinity of the receivers- However, an advantage 
of wider spacing is less sensitivity to noise. These 
competing concerns in combination with the wavelength 
of elastic or acoustic waves (or more accurately, the 
projection of the wavefield on the recording surface) , 

10 should be taken into consideration when designing the 
receiver group spacing. According to a presently 
preferred embodiment, the locally dense receivers are 
spaced about 1 meter apart. However spacing can in 
some situations be larger, for e.g. around 2 meters, or 

15 smaller, e.g. 0.5 meters. According to a preferred 
embodiment the inter-receiver spacing is around 0.25 
meters or less. As mentioned, the optimal spacing of 
the receivers depends upon the wavelength of interest. 
There should be at least two. receivers within the 

20 projection of the shortest wavelength of interest on 
the recording surface. According to a preferred 
embodiment the receivers are spaced apart by a distance 
approximately equal to or less than one-fifth of the 
shortest wavelength of interest, 

25 In the model based on the geometry shown in 

Figures 3a and 3b, the measurements from the receiver 
groups were used to obtain horizontal derivatives of 
the wavefield at each location. In this example, we 
will test our wavefield separation techniques on 

3 0 divergence. We expect that curl should display similar 
results. Finally, all sections shown in this example 
have been plotted using the same scale for amplitudes 
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so that amplitudes can be compared directly to each 
other both between traces and different sections. 

Figures 4a and 4b show the divergence as 
calculated directly by the reflectivity code. Figure 4a 
5 shows a section containing both up and down going 
waves. Note that as we approach zero offset the 
divergence vanishes. This is because the divergence is 
zero for vertically propagating plane P-waves. Figure 
4b shows the divergence of the up-going waves only. 

10 This is the desired P-wave section that we wish to 

obtain using our approach for wavefield separation and 
will therefore serve as our reference solution. Notice 
that the amplitude variation of different events in the 
section are quite different compared to those in the 

15 section in Figure 4a. Also notice the absence of some 
events due to mode conversions from S- to P-waves at 
the free surface. Some numerical artifacts in the 
numerical solution are also visible (e.g., the flat 
event before the first arrivals) . 

20 Traditionally, 3C data have been interpreted 

by assuming that waves propagate vertically near the 
receivers (steep gradients in material properties are 
assumed in the near-surface region) . Hence, P-waves 
show up on the vertical V3 component whereas S-waves 

25 appear on the horizontal vi and v 2 components. Figures 
5a-d illustrate sections obtained from V3 measurements 
only. To be able to compare them with divergence we 
have applied a time derivative to the v 3 measurements 
and scaled them with the P-velocity. Figure 5a shows v 3 

30 divided by a factor of two, according to a conventional 
technique. This exactly corresponds to the up-going P- 
waves for normal incidence (equation (36)). Figure 5b 
shows the difference between this section and the 
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reference solution. We see that the correspondence to 
P-waves rapidly breaks down away from normal incidence, 
where S-waves and mode conversions significantly 
contaminate the result. 
5 Part of the problem is of course that the V3 

measurements contain both up- and down going waves . 
According to the invention, these can be removed using 
a preferred wavefield separation technique as described 
herein. As can be seen from equations (31) and (35), 

10 this requires knowledge about P- and S-velocities in 

the near surface as well as convolving spatial filters 
with the measured vi and V2 components and adding them 
to the v 3 measurements. Figure 5c shows such a 
section. Figure 5d shows the difference between this 

15 solution and the reference solution. Even though the 
result has improved somewhat compared to just taking 
the raw v 3 measurements, the result can be improved 
even further . 

Figures 6a-d show the results of using a 

2 0 wavefield separation technique according to a preferred 
embodiment of the invention. In Figure 6a we have used 
the zeroth-order Taylor approximation (first term in 
equation (45)). This combines first-order spatial 
derivatives of vi and V2 with a time derivative of V3 . 

25 Figure 6b shows the difference between this solution 
and the reference solution. Although, the numerical 
noise that is present in the numerical solution has 
been amplified somewhat by the spatial derivatives, the 
result is now much closer to the reference solution. 
.30 In Figure 6c we have used the first-order 

Taylor approximation (both terms in equation (45)). 
This combines first-order spatial derivatives of vi and 
v 2 with a time derivative and second-order spatial 
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derivatives of v 3 . Figure 6d shows the difference 
between this solution and the reference solution. 
Again, the solution has improved even further without 
increasing the noise level compared to the zeroth-order 
5 appr oxima t i on . 

By comparing all the difference sections in 
Figures 5 and 6, it is clear that wavefield separation 
using the techniques according to the invention yield 
much better results than the traditional P-wave 

10 interpretation techniques. This is particularly the 
case for the window between 0.2 and 0.5 s and 0 m to 
150 m offset in the sections which contain events with 
very realistic incidence angles at moderate to low 
angles with respect to normal incidence (velocity 

15 gradients in the near surface causes energy to be 

incident fairly close to normal incidence) . In this 
example, we estimate that the approximate filters 
improve the results particularly for incidence angles 
up to 30° from normal incidence (this result is 

20 strongly dependent on material properties). However, it 
should be noted that only the theoretical expression 
for the divergence up/down separation filter in 
equation (40) is exact for all wave types. Therefore, 
longer filter approximations will tend to deal better 

25 with horizontally propagating and evanescent wave 
modes . 

According to the invention a new approach for 
P/S separation of land surface seismic data and 
removing the effects of the free surface is provided. 
3 0 By converting vertical derivatives to horizontal 
derivatives through the use of the free surface 
condition, the methodology can be used even when 
measurements are taken only at the free surface. 
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Therefore by making locally dense measurements of the 
wavefield, all spatial derivatives needed to compute 
divergence and curl of the wavef ield at the free 
surface can be obtained. These in turn correspond to P- 
5 and S-waves in isotropic media. 

The effects of the free surface can be 
removed through an up/down separation step as described 
herein. The filter for P-waves depends on both P- and 
S-velocity at the receivers, whereas the S-wave filters 

10 only depend on the S-velocity. Approximations to the 
filters were derived using Taylor approximations and 
tested on synthetic data. 

Filters have been derived for up/down and P/S 
separation using plane wave expansions in Dankbaar, J. 

15 W. M., 1985, Separation of P- and S-waves: Geophys. 
Prosp., 33, 970 - 986. These expressions can be 
compared to our expressions for full filters. A 
principal difference between the work by Dankbaar and 
those described herein, is that by deploying dense 

20 configurations of 3C geophones, P/S and up/down 

separation in 3D for each recording station can be done 
separately. The preferred up/down separation step makes 
use of approximations to spatial filters to obtain 
operators that are consistent with the number of 

25 geophones in the recording station. This approach is 
more robust since statics and near surface properties 
should be consistent within each recording station. 

The present invention can also be implemented 
in several steps where for instance the P/S separation 

3 0 is performed in 3D. The up/down separation can be 

performed in 3D or in 2D using an implicit filter if 
data are acquired along a 2D line for instance. 
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For 3C acquisition of land surface seismic 
data it is common practice simply to interpret the 
vertical component as the P-section and the horizontal 
components as SV- and SH-sections. This "traditional" 
5 P/S interpretation is exact for vertical arrivals. 

However, as energy is incident at non-normal incidence 
angles, this approximation tends to break down, both 
because the different waves appear on all components, 
but also because reflection coefficients differ from 

10 unity and mode-conversions occur at the free-surface. 
By comparing the "traditional" P-sections to the new 
methodology using synthetic data, we find a significant 
improvement in obtaining accurate amplitude and phases 
of arrivals for non-normal incidence angles . By simply 

15 using a zeroth-order Taylor approximation, we obtained 
sufficiently accurate results up to incidence angles of 
up to around 30° away from normal incidence (this 
result depends on material properties in the example) . 
Note that a zeroth-order Taylor approximation only 

2 0 involves first-order derivatives in time and space 

(along the free surface) . Note that the approximate 
expression for divergence consists of two terms. The 
first term corresponds to the divergence in the 
presence of the free surface scaled by a material 
25 constant. The second term is a time derivative of v 3 

scaled by a material constant. Hence, a correction is 
added to the "traditional" P-interpretation through the 
first term which improves the accuracy for incidence 
angles outside normal incidence. 

3 0 Figure 7 is a schematic illustration of a 

seismic data acquisition and processing system, 
according a preferred embodiment of the invention. 
Seismic sources 150, 152, and 154 are depicted which 
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impart vibrations into the earth at its surface 166. 
The vibrations imparted onto the surface 166 travel 
through the earth; this is schematically depicted in 
Figure 7 as arrows 170. The vibrations reflect off of 
5 certain subterranean surfaces, here depicted as surface 
162 and surface 164, and eventually reach and are 
detected by receiver groups 156, 158, and 160. 
Additionally, vibrations in the form of ground roll, 
shown as arrow 128, travel close the surface and can be 
10 recorded by the receiver groups. Each receiver group 
comprises a number of receivers such as in the 
arrangements depicted in Figure 3b. 



embodiment, the spacing of between the receivers within 
15 a single receiver group is substantially less than the 
spacing between the receiver groups. Schematically, 
this is shown in Figure 7 by the size 122 of receiver 
group 160 is substantially smaller than the distance 
120 between group 160 and an adjacent group 158. In 
20 the example shown in Figure 3a and 3b, the distance 120 
is 25 meters and the size 122 of the receiver groups is 
0 . 5 meters . 



receivers in groups 156, 158, and 160 convert the 
25 vibrations into electrical signals and transmit these 
signals to a central recording unit 170, usually 
located at the local field site. Preferably, the data 
is not group formed, but data from each geophone 
component are recorded. The central recording unit 
30 typically has data processing capability such that it 
can perform a cross-correlation with the source signal 
thereby producing a signal having the recorded 
vibrations compressed into relatively narrow wavelets 



Importantly, according to the preferred 



Referring again to Figure 7, each of the 
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or pulses. In addition, central recording units may 
provide other processing which may be desirable for a 
particular application. Once the central processing 
unit 170 performs the correlation and other desired 
5 processing, it typically stores the data in the form of 
time-domain traces on a magnetic tape. The data, in 
the form of magnetic tape is later sent for processing 
and analysis to a seismic data processing center, 
typically located in some other geographical location. 

10 The data transfer from the central recording unit 170 
in Figure 7 is depicted as arrow 176 to a data 
processor 180. Data processor 180 can be used to 
perform processing as described in steps 346 through 
352 as shown in Figure 8. 

15 Figure 8 shows steps involved in the 

filtering technique, according to preferred embodiments 
of the invention. In step 340 the seismic survey is 
designed. According to the present invention, the 
survey is preferably designed such that the receivers, 

2 0 for example geophones or hyrdophones, are positioned in 
locally dense arrangements as shown and described above 
with respect to Figure 3a, 3b, and 7. 



of wavefield separation is to be carried out. As 
25 described herein, analytical expression can be used for 
up/down separation, p/s separation or both. In the 
case of up/down separation, as described above, one 
would preferably use equations (29), (30), (31), (34) 
and (35) . In the case of p/s separation, as described 
30 above, one would preferably use equations (17) through 
(20) . In situations where both up/down and p/s 
separation is desired, one would preferably use 
equations (36) through (44) . 



In step 342 a choice is made as to which type 
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In step 346, an appropriate filter is 
designed based on the desired type of wavefield 
separation (and, accordingly, the equations selected) , 
and based on group recording geometry as designed in 
5 step 340. In creating an appropriate filter one should 
preferably choose a suitable approximation to fit the 
recording geometry. For example, a Taylor series 
approximation is described above, in equations (45)- 
(49). 

10 In step 344, the seismic data is measured 

transmitted, recorded and stored for each individual 
receiver in the survey. 

In step 350, a calculation is made of all 
partial derivatives in time and space needed for the 

15 filter approximations in step 346. 

In step 348, near surface material properties 
are estimated. This can be performed according to a 
number of methods. The presently preferred method is 
using techniques described in UK Patent Application 

20 entitled "System and Method for Estimating Seismic 
Material Properties" (UK Patent Application No. 
0003410.8) filed on 15 February 2000. However, one 
could also interpret near surface properties from the 
surface geology, or alternatively, a shallow refraction 

25 survey, or other suitable means could be used to 
estimate the near surface material properties. 

In step 352 the seismic data is filtered as 
determined in step 346. This is done by combining the 
estimated near surface material properties estimated in 

30 step 348, the seismic data recorded in step 344, and 

the partial derivatives of the wavefield calculated in 
step 350. The result is data 356 which is separated 
as desired (i.e. up/down, p/s or both). 
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While preferred embodiments of the invention 
have been described, the descriptions are merely 
illustrative and are not intended to limit the present 
invention. For example, while the preferred 
5 embodiments of the invention have been described 

primarily for use on the land surface, the invention is 
also applicable to receivers placed on and below the 
ocean floor. In the case of ocean bottom receivers, it 
is preferable to use stress conditions relevant for 

10 fluid-solid boundaries rather than the free surface 
condition. Additionally, the present invention is 
applicable to seismic measurements made in a borehole, 
known as borehole seismics. Although the examples 
described assume an essentially isotropic medium in the 

15 near surface region, the invention is also applicable 
to anisotropic media. In the case of anisotropic 
media, one may wish to increase the number geophones 
per group. 
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CLAIMS 



10 



15 



20 



25 



What is claimed is: 

1. A method of creating a filter for use with 
locally dense seismic data comprising the steps of: 
obtaining survey geometry characteristics 
from a locally dense seismic survey designed to 
record characteristics of an elastic or acoustic 
wavefield, the survey comprising a plurality of 
groups of receivers, and each group comprising at 
least three receivers densely spaced from each 
other; 

designing a filter which uses spatial 
derivatives of the wavefield of order between or 
including one and the maximum order of spatial 
derivatives of the wavefield that can be estimated 
within a group, such that the filter when 
combining data from within a single group, 
separates components of some or all of the 
wavefield arriving at the single group. 

2. A method according to claim 1 wherein 
the filter is designed to separate up/down going 
components of some or all of the wavefield. 

3. A method according to claim 1 wherein 
the filter is designed to separate p/s components of 
some or all of the wavefield. 

4. A method according to any of claims 1-3 
wherein the filter is designed to separate p/s and 
up/down components of some or all of the wavefield. 
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5 . A method according to any of the 
preceding claims wherein each of the densely spaced 
receivers in the group are spaced apart such that 
5 statics in the portion of the wavefield of interest are 
substantially constant. 



6 . A method according to any of the 
preceding claims wherein each of the densely spaced 
10 receivers in the group are spaced apart by about 2 
meters or less. 



7 . A method according to any of the 
preceding claims wherein each of the densely spaced 

15 receivers in the group are spaced apart by a distance 
of about one fifth the shortest wavelength of interest 
or less. 

8 . A method according to any of the 

20 preceding claims wherein the locally dense seismic 

survey comprises generating elastic or acoustic waves, 
and the receivers in a group span less than about the 
smallest wavelength of said elastic or acoustic waves. 

25 9. A method according to claim 1 further 

comprising the step of calculating partial derivatives 
of the wavefield. 

10. A method according to claim 9 wherein 
30 the step of calculating partial derivatives includes 
using a taylor series expansion as an approximation. 
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11. A method according to claim 9 further 
comprising the step of estimating near-surface material 
properties . 



5 12 . A method according to claim 11 further 

comprising the step of filtering 'the seismic data by 
combining the estimated near surface material 
properties, the seismic data, and the calculated 
partial derivatives. 

10 

13 . A method according to claim 1 wherein a 
portion of the wave field that is separated is the 
pressure wavefield. 



15 14 . A method according to claim 13 wherein 

the filter separates surface waves and/or air wave 
induced ground motion from the seismic data. 

15 . A method according to claim 1 wherein 
20 the locally dense seismic survey is performed on land. 

16. A method according to claim 15 wherein 
the step of creating the filter comprises using the 
free surface condition to convert vertical derivatives 

2 5 of the wavefield to horizontal derivatives of the 
wavefield 



17 . A method according to claim 1 wherein 
the seismic survey is performed primarily for 
30 hydrocarbon reservoir exploration, evaluation or 
characterisation . 
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18. A method according to any of the 

preceding claims wherein the near surface velocity is 
essentially isotropic. 

5 19 - A method according to any of the 

preceding claims wherein the near surface velocity is 
anisotropic . 



WO 01/53854 PCT/GB01/00186 



1/15 




Figure 1a 




Figure 1 c 



WO 01/53854 PCT/GB01/00186 



2/15 



5( : ~— l — I r— — r T 




,| ■ , 1 ^ 1 ■■! » 1 

0 O.flS H.1 0,t5 0,2 L\2B 4*3 

Itoiteatttai viavtntiirtjor 



Figure 2 



WO 01/53854 PCT/GB01/00186 



3/15 



204 




T.TTTTTTT 



450 m 



212 



Vj = 1500 m/s 

r>- f - H 500 m/s C; 



Vpia.lSOOtn/s"' 
900 wis 



Figure 3a 



150 re 



210. 



250 



256 



0.5m » 



252 



-254 



T 



260 



Figure 3b 



WO 01/53854 PCT/GB01/00186 



4/15 




Figure 4a 



WO 01/53854 PCT/GB01/00186 



# 



5/15 




Figure 4b 



WO 01/53854 

# 



6/15 



PCT/GBOl/00186 



v^tt. Surface fMsaaiueTOfiiitft 




Figure 5a (prior art) 



WO 01/53854 PCT/GB01/00186 



7/15 




Wtfer&f»a :ttfw?fin rdkenoo v & 




eflscl ^mj 



Figure 5b (prior art) 



WO 01/53854 



8/15 



PCT/CB01/00186 




Figure 5c 



WO 01/53854 PCT/GB01/00186 




9/15 




WO 01/53854 PCT/GB01/00186 



10/15 




Figure 6a 



WO 01/53854 PCT/GB01/00186 



11/15 



<v . v) u Cutanea between reler ence dud am, outer appjte 




Figure 6b 



WO 01/53854 PCT/GB01/00186 



12/15 





Figure 6c 



WO 01/53854 PCT/GB01/00186 



13/15 





Figure 6d 



WO 01/53854 PCT/GB01/00186 



14/15 






WO 01/53854 



PCT/GB01/00186 




15/15 



Design survey and 
locally dense group 
recording geometry 

340 



T 



Choose up/down and/or 
P/S separation 

342 



Record data 
344 



Estimate near- 
surface properties 

348 



Filter design 
346 



Calculate partial 
derivatives of 
wavefield 

350 



Filter data 



352 



I 



Up/down and/or P/S 
separated data 

356 



Figure 8 



INTERNATIONAL SEARCH REPORT 



A. CLASSIFICATION OF SUBJECT 

IPC 7 G01V1/36 



GtllV 



I onal Application No 

PCT/^fcl/00186 



1V1/28 



According to International Patent Classification (IPC) or to both national classification and IPC 



B. FIELDS SEARCHED 



Minimum documentation searched (classification system followed by classification symbols) 

IPC 7 G01V G10K 



Documentation searched other than minimum documentation to tho extent thai such documents arc included in the fields searched 



Electronic data base consulted during the International search (name of data base and, where practical, search terms used) 

EPO-Internal , INSPEC, COMPENDEX, WPI Data 



C. DOCUMENTS CONSIDERED TO BE RELEVANT 



Category " Citation of document, with indication, where appropriate, of the relevant passages 



GB 2 333 364 A (GECO AS) 

21 July 1999 (1999-07-21) 

cited in the application 

page 7, line 10 -page 13, line 21 

EP 0 327 758 A (CONOCO INC) 
16 August 1989 (1989-08-16) 
abstract 

AMUNDSEN L: "Wavenumber-based filtering 

of marine point-source data" 

GEOPHYSICS, SEPT. 1993, USA, 

vol. 58, no. 9, pages 1335-1348, 

XP002166869 

ISSN: 0016-8033 

abstract 

page 1338, column 2, line 26 -page 1340, 
column 1, line 33 

-/-- 



Relevant to dalm No. 



Further documents are listed In the continuation of box C. 



[X I Patent family members are listed In annex. 



° Special categories of cited documents : 

•A* document defining the general state of the art which is not 
considered to be of particular relevance 

"E" earlier document but published on or after the International 
filing date 

V document which may throw doubts on priority claim(s) or 
which is cited to establish the publication date of another 
citation or other special reason (as specified) 

"O" document referring to an oral disclosure, use, exhibition or 
other means 

*P* document published prior to the international filing date but 
later than the priority date claimed 



•T" later document published after the international filing date 
or priority date and not in conflict with the application but 
cited to understand the principle or theory underlying the 
invention 

"X* document of particular relevance; the claimed invention 
cannot be considered novel or cannot be considered to 
Involve an inventive step when the document is taken alone 

"Y" document of particular relevance; the claimed invention 
cannot be considered to involve an inventive step when the 
document is combined with one or more other such docu- 
ments, such combination being obvious to a person skilled 
in the art. 

"&' document member of the same patent family 



Date of the actual completion of the International search 



25 May 2001 



Date of mailing of the International search report 



08/06/2001 



Name and mailing address of the ISA 

European Patent Office, P.B. 5818 Patentlaan 2 
NL-2280 HV Rijswijk 
Tel. (+31-70) 340-2040, Tx. 31 651 epo nl, 
Fax: (+31-70) 340-3016 . 



Authorized officer 



Lome, B 



Foim PCT/ISA/210 (second shoot) {July 1992) 



page 1 of 2 



^Continuation) DOCUMENTS COI 



'RED TO BE RELEVANT 



I t^^Appllcation No 

PC TjJbl /00186 



Category 0 Citation of document, with indicallon.where appropriate, of the relevant passages 



Relevant to claim No. 



ROBERTSSON J 0 A ET AL: "Viscoelastlc 

finite-difference modeling" 

GEOPHYSICS, SEPT. 1994, USA, 

vol. 59, no. 9, pages 1444-1456, 

XP002166870 

ISSN: 0016-8033 

abstract 

page 1448, column 2, line 26 -page 1449, 
column 2, line 11 

AMUNDSEN LASSE ET AL: "Decomposition of 

multi component sea-floor data into upgoing 

and downgoing P- and S-waves" 

GEOPHYSICS, US, TULSA, OK, 

vol. 60, no. 2, March 1995 (1995-03), 

pages 563-572, XP002130603 

abstract 

page 565, column 1 -column 2 

ZHOU B ET AL: "Numerical seismogram 
computations for Inhomogeneous media using 
a short, variable length convolutional 
differentiator" 

GEOPHYSICAL PROSPECTING, AUG. 1993, 
NETHERLANDS, 

vol. 41, no. 6, pages 751-766, 
XP000994570 
ISSN: 0016-8025 
page 753 -page 754 



Form PCT/ISA/210 (continuation of cooond enact) (July 1092) 



page 2 of 2 



Patent document 
cited in search report 



an on patent family members 


tlong^jllcatlon No 

PCT/G^h)0186 


Publication 
date 


Patent family 
member(s) 


Publication 
date 



GB 2333364 



21-07-1999 



FR 
NO 
US 



2773618 A 
990153 A 
6101448 A 



16-07-1999 
16-07-1999 
08-08-2000 



EP 0327758 


A 


16-08-1989 


US 


4794573 A 


27-12-1988 








CN 


1035002 A,B 


23-08-1989 








OP 


1248087 A 


03-10-1989 








NO 


884738 A 


14-08-1989 



I 

i 



Form PCT/ISA/210 (patent family annex) (July 1992) 



